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We study the influence of difTerent edge types on the electronic density of states of graphene 
nanostructures. To this end we develop an exact expansion for the single particle Green's function 
of ballistic graphene structures in terms of multiple reflections from the system boundary, that 
allows for a natural treatment of edge effects. We first apply this formalism to calculate the average 
density of states of graphene billiards. While the leading term in the corresponding Weyl expansion 
is proportional to the billiard area, we find that the contribution that usually scales with the total 
length of the system boundary differs signiflcantly from what one flnds in semiconductor-based, 
Schrodinger type billiards: The latter term vanishes for armchair and infinite mass edges and is 
proportional to the zigzag edge length, highlighting the prominent role of zigzag edges in graphene. 
We then compute analytical expressions for the density of states oscillations and energy levels within 
a trajectory based semiclassical approach. We derive a Dirac version of Gutzwiller's trace formula 
for classically chaotic graphene billiards and further obtain semiclassical trace formulae for the 
density of states oscillations in regular graphene cavities. We find that edge dependent interference 
of pseudospins in graphene crucially affects the quantum spectrum. 

PACS numbers: 73.22.Pr, 73.22.Dj, 73.20.At, 03.65.Sq 



I. INTRODUCTION 
A. Graphene-based nanostructures 

Triggered by the experimental discovery of massless 
Dirac quasiparticlesi*^, graphene has become one of the 
most intensively studied materials of the last decade (for 
reviews on physical properties see Refs. 0-0) • 

Subsequently, graphene-based nanostructures have 
been the focus of an immense experimental activity, in- 
cluding graphene nanoribbons^"— , quantum dots^^— 



.,17.18 



rais- 



Aharonov-Bohm ring s i and antidot array: 
ing the issue of confining massless Dirac electrons. On 
the theoretical side, several studies have also focused on 
graphene nanostructures: Graphene nanoribbons have 
been studied first using a lattice modeU^i^. The wave- 
functions and energy spectra of graphene nanoribbons 
have been derived by Brey and Fertig^i for armchair and 
zigzag type edges, and by Tworzydio and coworkers^ for 
the case of infinite mass edges. The spectral and trans- 
port properties of Dirac electrons confined in graphene 
quantum dots have been investigated analyticallj«2^"— 
and by numerical means^^^— . Also energy spectrum and 
conductance of Aharonov-Bohm rings have been the fo- 
cus of several publications^"— as well as superlattice ef- 
fects in graphene antidot lattices^i^ and the density of 
states of nanoribbon-superconductor junctions^. 

One upshot of these studies is the understanding that 
the confinement of charge carriers in graphene affects the 
coherent electron and hole dynamics considerably. In 
conventional two-dimensional electron systems (2DES) 
such as low-dimensional semiconductor structures, the 



charge carriers can be confined, e.g. by the application of 
top or side gate voltages, and the quasiparticle transport 
does not depend on the minute details of the resulting ef- 
fective potential. In contrast, in graphene, electrostatic 
potentials do not necessarily confine charge carriers as 
the Dirac spectrum does not have a gap*^. Thus the con- 
fined electrons or holes in graphene nanostructures or 
flakes are expected to scatter from the very ends of the 
terminated graphene lattice, and the internal degrees of 
freedom (such as spin or pseudospin) of the quasiparticles 
before and after the scattering are considerably affected 
by the atomic level details of the edges. This mixing of 
internal (pseudo)spin with orbital degrees of freedom of 
charge carriers at the boundary leads to richer boundary 
conditions than for the conventional 2DES^"— . These 
boundary conditions in turn affect the spectral and trans- 
port properties. However, experimental control and ma- 
nipulation of edges at an atomistic level is far from be- 
ing achieved. Thus a full theoretical description is de- 
sirable. However, the edge disorder differs from usual 
(weak) bulk disorder in that weak coupling perturbation 
theories cannot treat edges. Therefore this paper is ded- 
icated to develop a formalism that includes the effects of 
edges non-perturbatively, and to subsequently apply this 
formalism to study edge effects on the spectral density of 
states of graphene nanostructures. 



B. Scope of this work 

Cutting a finite piece of graphene out of the bulk 
will generally lead to disordered boundaries with local 
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properties depending on the respective orientation of an 
edge segment with respect to the crystallographic axes. 
The accurate calculation of the eigenenergies these fi- 
nite graphene systems usually requires numerical quan- 
tum mechanical approaches. However, it appears diffi- 
cult to systematically resolve edge phenomena from other 
quantum effects or to unravel generic features of graphene 
nanostructurcs using numerical simulations. Here we fol- 
low a complementary strategy: We adapt the multiple 
reflection expansio n^^''*" , i. e. a representation of the 
Green's function in terms of the number of reflections 
from the system boundaries, to the case of graphene. We 
thus incorporate edge effects (due to armchair, zigzag 
and infinite mass type and combinations of such edge seg- 
ments) in a direct and transparent way. We next derive 
a semiclassical approximation for the Green's function, 
assuming the Fermi wavelength is much smaller that the 
typical system size i, i. e. L ^ l/fce. On the other hand, 
the Dirac equation that we use is valid for Fermi wave- 
lengths that are large compared to the lattice constant 
aw 2. 46 A, i.e. ifl/fce^a. For mesoscopic systems 
with L ^ a, the semiclassical approximation can thus 
be well fulfilled in the linear dispersion regime, in which 
quasiparticlc dynamics is governed by the effective Dirac 
equation. The resulting Green's function then can be 
used to calculate the density of states (DOS) or the con- 
ductance, and their correlators. 

In this work we consider the density of states. Wc 
focus on gross structures and spectral densities arising 
from moderate smearing of the level density and on the 
calculation of DOS oscillations and individual levels sep- 
arately. To this end we decompose the DOS into an 
average part and the remaining oscillatory contribution. 
The average spectral density, approximated by the so- 
called Weyl expansion^Sii^i^ valid in the semiclassical 
limit, is a fundamental quantity of a cavity. It incorpo- 
rates various geometrical and quantum features, includ- 
ing edge effects. For billiards with spin-orbit interaction, 
the smooth part of the engery spectrum has been studied 
in Ref.Sl- The oscillatory part of the DOS is computed 
by invoking a semiclassical approximation, leading to so- 
called semiclassical trace formulae, i. e. sums over co- 
herent amplitudes associated to classical periodic orbits. 
For graphene cavities with shapes giving rise to regular or 
chaotic classical dynamics we derive trace formulae anal- 
ogous to those known (Berry- Tabor— and Gutzwiller— 
formula, respectively) for the corresponding Schrodinger 
billiards, i. e. billiard systems based on the Schrodinger 
equation with Dirichlet boundary conditions. For two 
representative regular shapes, we compute the DOS os- 
cillations and the semiclassical energy levels explicitly. 
The effects of both, the underlying effective Dirac equa- 
tion (for graphene close to the Dirac point) and reflec- 
tions at different kinds of edges, is incorporated by a 
pseudospin propagator associated with each orbit, mul- 
tiplying the usual semiclassical amplitude. Semiclassi- 
cal trace formulae involving the electron spin dynamics 
have been earlier considered for the massive Dirac equa- 



tion by Bolte and Keppeler— and for bulk graphene by 
Carmier and Ullmc4^. Related trace formulae appear 
also in trajectory-based treatments of electronic systems 
with spin-orbit interaction^"—. We note that semiclas- 
sical methods have also been used to study graphene in 
magnetic fields^"—. 

Following the concepts outlined above wc address edge 
effects on the electronic spectra of closed graphene cavi- 
ties and quantum transport through open graphene sys- 
tems in two consecutive papers. In the present paper 
we first derive the single-particle Green's function and 
its semiclassical approximation for graphene cavities and 
calculate the density of states. In subsequent work2& 
we will consider quantities based on products of single- 
particle Green's functions. They include the transport 
quantities such as the conductance as well as the spec- 
tral two-point correlator and its dual the spectral form 
factor, as a tool to study spectral statistics. The semi- 
classical treatment of observables based on products of 
Green's functions requires additional techniques which 
builds the conceptual basis of the second paper—. 

The present paper is organized as follows: After in- 
troducing below the effective Hamiltonian and (matrix) 
boundary conditions for the different edge types, we de- 
rive in Sec. |lT]the multiple reflection expansion (MRE) 
for the Green's function of a ballistic graphene structure. 
With this expansion as a starting point we then compute 
in Sec. IIIII the first two terms in the Weyl expansion 
for the smooth part of the DOS of graphene billiards, 
particularly focusing on contributions from the bound- 
ary. We compare our analytical theory with numerical 
quantum simulations for various graphene billiards with 
different edge structures. In Sec. IIVI we turn to the os- 
cillatory part of the DOS . To this end we first obtain a 
general semiclassical approximation to the MRE for the 
graphene Green's function in terms of sums over classical 
trajectories in IIV Al Subsequently we focus on the DOS 
oscillations in graphene billiards with regular classical dy- 
namics in IIVBI Wc give semiclassical trace formulae for 
two exemplary geometries, namely disks and rectangles, 
and discuss the effects of the graphene edges. Finally we 
extend Gutzwiller's trace formula for the oscillatory part 
of the DOS to graphene cavities with chaotic classical 
dynamics in IIV Gl We conclude in Sec. |V] and gather 
further technical material in the appendices. 

C. Hamiltonian and boundary conditions 

Neglecting the conventional spin degree of freedom, the 
effective Hamiltonian that describes electron and hole dy- 
namics in graphene close to half filling is^ 

H = VpT;: ® CTx + VpTo <Sl (Ty Py , (1) 

where Vp is graphene's Fermi velocity. The {at} de- 
note Pauli matrices in sublattice pseudospin space and 
Pauli matrices in valley-spin space are repescnted by 
{Ti}, while (To a-nd tq are unit matrices acting on the 
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corresponding spin space. In the following, we usually 
omit the latter. The Hamiltonian ([1]) acts on spinors 
[ipA, '0A' V'b] where A/B stands for the sublattice in- 
dex and the primed and unprimed entries correspond 
to the two valleys. We find it convenient to transform 
Eq. ([T]) to the valley isotropic formal using the unitary 
transformation 

1 i 
^ i^{Ta + T,) ® + -(t^ - T,)ay . (2) 

The transformed Hamiltonian is 

H = U^HU ^VpTo®(T-p (3) 

and acts on spinors [V'a, V'Bi ""05) V'a]- 

We consider a graphene flake in which electron and 
hole dynamics is confined to an area V. The boundary 
condition on the spinors at a point a on the boundary 
dV is expressed as Paip\a = 0, where Pa is a 4 x 4 
projection matri x"^^'^^ . Throughout this paper we reserve 
bold Greek letters for boundary points and bold Roman 
letters for points in the bulk of the flake. For the most 
common boundaries, i.e. zigzag (zz), armchair (ac) and 
infinite mass (im), the boundary matrices are given hy^ 

P„ = i(l-iy-x®J7-cr) (4) 

where the vectors v and t] are summarized in Tab. ID 
K = 4ir/3a is the distance of the Dirac points from the 
F-point of the reciprocal space, Xa = a ■ x and ta is the 
direction of the tangent to dV at a. For zigzag edges 
the sign in r/ is determined by the sublattice of which 
the zigzag edge consists. For an A-edge the upper sign 
is valid and for a i?-edge the lower sign. That means 
the orientation of the edge effectively determines rj. For 
armchair edges, the upper sign is valid when the order of 
the atoms within each dimer is A-B along the direction 
of ta, and the lower sign is valid for B-A ordering. For 
infinite mass edges, the sign depends only on the sign of 
the infinite mass. The upper sign is valid for the mass 
going to +0O outside of V and the lower for the mass 
going to — cx). 

We note that for a model that includes next near- 
est neighbour hopping (nnn), the boundary conditions 
need to be modified to include differential operations on 
the spinor. Nevertheless, as we shall show in App.[Bl it 
is possible to modify our formalism to account for nnn 
hopping approximately by keeping only nearest neighbor 
hoppings, but modifying the boundary conditions intro- 
ducing an edge potential. 

D. Single particle density of states 

The single particle DOS for a closed system is defined 
as^ 

p{k^)=J2HkE~K) ■ (5) 

n 
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im 
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+ cos{2Kxo,)y 
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TABLE I. The vectors u and t/ for zigzag (zz), armchair (ac) 
and infinite mass (im) type boundaries. 

Here n labels the eigenenergics i?„ = hvpkn, and we de- 
fine E = hvpkE- In our derivation below we use the rela- 
tion between the DOS and the retarded Green's function 
of a system, 

pikp) = -;^3m j dx Tr [G{x, x)] , (6) 

V 

where the Green's function G fulfills 

{E + iT]-H)G{x,x')^hvpS{x-x'), (7) 

with the Hamiltonian H acting on the first argument 
of G. For a mesoscopic graphene flake the mean level 
spacing Ak, which is given by the inverse area of the 
system, is typically of the order 10~* 1/a or smaller. This 
means that p is in principle a rapidly oscillating function 
of kp. However, one can decompose p into a smooth part 
p and an oscillating part posc in a well defined way^i^, 

P= P + Pose ■ (8) 

In this work, wc address both contributions to p and fo- 
cus on the particularities that arise due to the spinor 
character and the linear dispersion of quasiparticles in 
graphene. The smooth part p represents the density of 
states in the limit of strong level broadening. Technically, 
level broadening is achieved by adding a finite imaginary 
part to the Fermi energy or in other words considering 
a real self energy. This corresponds to an exponential 
damping of the Green's function and therefore only tra- 
jectories of short length, in the limiting case of 'zero- 
length', contribute. In Sec. IIIII we treat p in detail. On 
the other hand, posc is connected to (periodic) orbits of 
finite length, and in Sec. IIVI we use a semiclassical ap- 
proach to describe this part of the density of states. 

In the following, we derive an exact expression for 
the Green's function entering Eq. ^ and later also its 
asymptotic form in the semiclassical limit, valid for large 
system sizes. 



II. THE MULTIPLE REFLECTION EXPANSION 
FOR GRAPHENE 

In this chapter, we derive a formula for the exact 
Green's function of a graphene cavity. The Green's func- 
tion can then be used to obtain e. g. the spectral density 
of states or the conductance. In addition to Eq. G 
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FIG. 1. Schematic representation of a quantum path con- 
tributing to the Green's function G{x,x'). The black hues 
with arrows stand for free propagations described by Go , while 
each black disk represents a vertex of the form j(T„^ Pa ■ 

also obeys the boundary conditions PctG(a,x') ~ for 
any given point a on the boundary. 

We now parameterize the full Green's function as a 
sum of the free retarded Green's function Go of extended 
graphene and a boundary correction that is produced by 
a, yet unknown, Dirac-chargc layer fi: 

G{x,x') = Go{x,x')-Jdaf3Goix,l3)ian^ fj,{(3,x') . (9) 
dv 

Here = cr v for an arbitrary vector v, and np stands 
for the normal unit vector at the boundary point (3 point- 
ing into the interior of the system. The free Green's 
function is obtained by solving Eq. ([7]) with boundary 
conditions Go{x, x') as |a; — x'\ -H- oo. It is given by 

Ga{x,x') = hvp{x\{E - ny^lx') 

= -^{k^ - iV^-(7)H+{k^\x - x'\) , (10) 

where Hq denotes the zeroth order Hankcl function of 
the first kind. The free Dirac Green's function can be 
expressed in terms of the free Schrodinger Green's func- 
tion go as 

Go{x,x') = (fefi - iVa;-cT)go{x,x') . (11) 

The Schrodinger Green's function t/o is a solution to 

{kl + iij - pyfe)go{x, x') = S{x - x') . (12) 

The parametrization in Eq. ([9]) is singular in the limit 
x^a^: 

lim G{x,x') = Gn{oi.,x') — —ii{a,x') (13) 
-J dai3Goict,(3)ianp f-i{f3,x') . 

dV 

The source of this singular behavior is the logarithmic 
divergence of -ffo^(C) as ^ 0. For a detailed derivation 
of Eq. (in]) sec App.S Multiplying ^ with and 



invoking the boundary conditions, we obtain an inhomo- 
geneous integral equation for the charge layer /i. As a 
first step we assume that Paji = /i, so that wc get 

^l{a,x')=2P^Goia,x') (14) 
-2 J dai3PaGo{a,f3) ian^ fj,{(3,x') . 

dV 

Since P^ ~ Pa, the unique solution of Eq. p^ . obtained 
by iteration, automatically fuUfiUs Pq./^ = and thus 
is already a solution of the original integral equation for 
^. Substituting this solution into Eq. (jH), we obtain the 
following expansion for the exact Green's function of a 
graphene fiake with generic edges: 

oo 

G{x,x')^G^{x,x')+^Gn{x,x'). (15) 

where 

Gn{x, x') = (-2)^ j d(Jo,M ■ ■ ■ daa^daa, x (16) 
dv 

Go{x,CXN)i<Jn^,^Pa!^ ■ ■ • Gq ((12 , «! )io-„ci -^ai Gq (oil , a;') . 

Each term in this expansion can be viewed as a se- 
quence of free propagations connected at reflections at 
the boundary (sec Fig. [I}. Wc thus obtain the multiple re- 
flection expansion (MRE). In Eq. (|16p every reflection is 
represented by a boundary dependent projection Pa and 
by (T„^ , a reflection of the pseudospin across the normal 
axis given by ria ■ The integrals along the boundary can 
be interpreted as a 'summation' over all quantum paths 
leading from x' to x. In Fig.[l] we show schematically a 
typical term in the MRE using the example of a quantum 
path that includes three reflections at the boundary. To 
summarize at this stage, with Eqs. (|15[ I16p wc obtained 
a formalism that naturally relates the edge effects to any 
quantity that involves single particle Green's functions. 

III. THE SMOOTHED DENSITY OF STATES 
OF GRAPHENE BILLIARDS 

A. Weyl expansion 

In the following we are going to derive the leading order 
contributions to the smoothed density of states p. In 
usual Schrodinger billiards of linear system size L, as 
they are realized e. g. in 2DES in GaAs heterostructures, 
p can be expanded in powers of ksL with leading order 
(ksL)^, a constant term (kEpy^ and higher order terms 
{ksL)"^, {kEL)~'^ and so forth as 

P = Po + Pi + P2 + P3 . . . . (17) 

In the large k^L limit, p is dominated by the first term, 
which docs not depend on the shape of the system but 
only on its total area. This theorem goes back to Her- 
mann Weyl^ and therefore the series is known as the 
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dV 




FIG. 2. For the calculation of the one-reflection term in 
the expansion for p, we work in the plane approximation: 
For a given point ct on the boundary dV, we approximate 
the boundary locally by the tangent at a and introduce a 
local coordinate system with x and y along the tangential 
and normal direction respectively. 

Weyl expansion for the density of states. Each of the 
terms in Eq. ((T7| can be obtained from the MRE (|16p : po 
originates from the zero-reflection term (simply Gq) and 
therefore scales with the total area A of the system. The 
term pi is due boundary contributions, obtained within 
the so-called plane approximation (cf. Fig. [2]), leading to 
a scaling with the length of the boundary. The term p2 
stems from curvature and corner corrections to the plane 
approximation and so forth. In this work we focus on 
leading contributions po and pi. The smooth contribu- 
tions are of qualitatively different origin than the oscil- 
lating part of the DOS, treated in Sec. IIVI While the 
latter correspond to orbits for which the phases occur- 
ing in Eq. ([6]) are stationary, the smooth DOS is due to 
trajectories approaching 'zero-length' for which the am- 
plitudes diverge. We find that the linear term in the Weyl 
expansion for graphene po is similar to the usual 2DES 
case, but the term pi behaves strikingly different. 

B. Bulk term 

We begin with the zero-reflection term Go{x,x) in 
graphene. From Eq. ([TU]) we can directly read off 

TT[Gaix,x')] = -ik^H+{k^\x - x'\) . (18) 

Although Go diverges as x' ^ Xr^ its imaginary part is 
finite. We get 

3mTr[Go(a;,a;)] = -Ifc^l . (19) 

Since there is no x dependence left, the spatial integral 
in Eq. ([5]) gives just A = | V|, the area of the billiard, and 
we have 

PoikE) ^ -\k^\ . (20) 

As for Schrodingcr billiards, the bulk term is pro- 
portional to the total area of the system. The energy 



dependence of po is however different, since fc^^ scales lin- 
early with energy in graphene but has a square root de- 
pendence in the Schrodinger case. 

C. Boundary term 

1. Plane approximation 

As we show below, the boundary term pi depends on 
as well as on the boundary length of the system, in 
a manner distinctly different from that of Schrodinger 
billiards. In order to evaluate pi, we assume that the 
energy has a finite imaginary part ^. This smoothcns the 
DOS and makes Gq an exponentially decaying function 
of the distance between x and x'. We start from Eq. ([5]), 
omit the free propagation term that led to po, and obtain 
for the remaining contribution to the smooth DOS 

(5p=— J daaj dxTr[Go{x,a.)ian^Pi{a,x)] . 

(21) 

Here we replaced the boundary integration by a sum of 
integrations over boundary pieces 9Vi, where the bound- 
ary condition is constant for each i. Further pi{a,x) is 
defined via Eq. ([T4| with a £ dVi . Since Go is short 
ranged, the dominant contribution to the boundary inte- 
gral in Eq. (|2ip comes from configurations where x is near 
the boundary point a, and the integral in Eq. (|14|) is dom- 
inated by contributions where f3 is near a. Thus we ap- 
proximate the surface near a by a plane (cf. Fig. [2]). The 
corrections to this approximation are of order l/k^R, 
with the local radius of curvature R ^ L, thus of higher 
order in the Weyl expansion^. We now take advantage of 
the homogeneity of the approximate surface at a and use 
Fourier transformation along the direction of the tangent 
to the dVi at a, to get for Sp sa pi 

OO CJO 

Pi = l^raJ2 m\ jdy, J ^Tr [SG,{k, y,)] , (22) 

* -OO 

with 

5Gi{k, yi) = Go(fc, yi)i(Jrt^pi{k, y^) . (23) 

Here yi is the ordinate of the local coordinate system at 
OL (see Fig. (2) and 

M,(fc,2/,) = 2F,(fc)P„Go(fc,-y.), (24) 

T^{k) = [l + 2P„Go(fc,0)iay]-' , (25) 

with the Fourier transform defined as 

OO 

/(^,2/)= / §e'^-/(fc,j/). (26) 

— OO 

We pushed the upper limits of the ?;,-intcgration to infin- 
ity, which is valid when cxp[— CJm(fcE)_L] ^ 1. To obtain 
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Eq. (|22p , we further assumed that a. is away from the 
corners where the boundary condition changes. The cor- 
rections due to such points arc of order l/ksL smaUer 
than the boundary term. 

The free Green's function in mixed representation is 
given by 



Go(fc,yj) 



_g-a(fc)|j/i 

2a(fc) 



with 



a{k) = v/fc2 - fc|, 5ne[a(A:)] > . 



(27) 



(28) 



Next we focus on contributions to the boundary term 
from various types of edges. 



2. Zigzag edge 
For a zigzag edge (without nun hopping, sec Tab. |I| 



Pc. = (lTT-^®CT^)/2. 



(29) 



Then F,; is diagonal in valley space and we can invert the 
valley subblocks separately giving 

r,(fc) - [a{k) - (fca, - ik,Oy){l - P„)] . 

(30) 

We insert Fi(fc), Eq. into Eq. dMl) and take into ac- 
count that Pa is a projection matrix, i.e. = Pa, to 
obtain for the Dirac-charge density 

,a(fc) 

kl 



/x,(fc,2/,) = -2-^[a(fc)±A:T,]P„Go(fc,-y,)- (31) 



Substituting this expression into Eq. (j23p . we obtain 
8G,{k,y,) (32) 

^ [a{k) ± fcr,] Go(fc, y,) ia^ P„ Go(fc, -y,) . 

Then the trace is given by (note that yi > 0) 

2fc2 



Tr[<5G,(fc,y,)] = - 



a(k)ki, 



-2a{k)yi 



(33) 



Evaluating the y^-integral we get (note that the real part 
of a(fc) is positive) 

oo 

/f dk 
dy^ j — Tr [<5G,(fc, y,)] = K^^M^b) , (34) 



where 



1 i 



(35) 



and we have introduced a cut-off momentum fcmax ~ !/«• 
Such a cut-off is justified, since in real graphene the avail- 
able fc-space is not infinite owing to the lattice structure. 



We cannot calculate the precise numerical value for fcmax 
within our effective model. Using tight-binding calcu- 
lations we estimate fcmax ~ n/Sa^. The result ([34|) 
means that without nun hopping, zigzag edges lead to 
a DOS contribution that is strongly peaked at zero en- 
ergy. The origin of this contribution is indeed the exis- 
tence of zigzag edge states at zero energyi^ii^i^i^i^. To 
understand this connection we consider the prefactors in 
Eq. ([31]) and Eq. in the limit of fc^ 0; then we 
have 



a(fc) 

1.2 



a{k) ± kTz 



1,2 



[1 ± sgn(fc) Tz 



(36) 



For the upper sign, this expression is divergent in one 
valley for negative k (r = +1) and in the other valley 
for positive k {t = —1) as ks approaches zero. For the 
lower sign it is just vice versa. Thus we identify the zero- 
energy states that are localized at the zigzag graphene 
edge. In a single valley this causes a strong asymmetry 
in the spectrum and breaks the (effective) time reversal 
symmetry. Below we show that the zigzag edge states 
are the only contribution to the DOS that scales with 
the boundary length of the graphene flake. Armchair and 
infinite mass type edges do not contribute to the surface 
term. However for the zigzag edge states, the effect of 
nnn hopping is significant22i^i^. For a more realistic 
description of the their effects on the DOS, it is therefore 
necessary to consider nnn hopping for the boundary term 
at zigzag edges. In App.|B] we show that the boundary 
condition for zigzag edges is effectively modified due to 
nnn hopping resulting in a boundary matrix 

P„ = i (1 T ® - ^t'cTy ± t'Tz (g) a^) . (37) 

Here t' 1 is the ratio of the nnn hopping integral and 
the nearest neighbor hopping integral in the tight-binding 
formalism. The effect of this boundary condition is to 
modify Eq. to 

Kk, y') = Mk) ^}!:l''[t^^Z Po^ Go{k, -y') ■ 



[a{k)-t'k^]'^-k^' 



(38) 



Note that the Eqs. (|37l [38 |) turn into the expressions 

for t' = 0. Following the same line of calculation we 

find 



Tr[5G,(fc,y,)] 



2fc2 



a{k) (1 - i'2)fc^ + 2<'a(fc) 



-2a(fc)y. 



(39) 



and the corresponding contribution to the DOS is to lin- 
ear order in t' 



7 J ^Tr[5G,(fc,2/,)] « 

-oo 



Here 



Of (fcij) = - arctan(fcE/^) + \ 

TT Z 



(40) 



(41) 
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is a smooth approximation to the Heaviside step function. 

According to Eq. PO)) . the fce-dependence of the zigzag 
contribution to the DOS is quahtatively altered by the 
inchision of nnn hopping. It is strongly asymmetric due 
to the broken electron-hole symmetry^. Also the peak at 
zero = has disappeared, because the edge states are 
not degenerate anymore but exhibit a linear dispersion 
^cdgc _ ^^2t' as derived in App. |B] Note that in tight- 
binding, there is still a van Hove singularity in the DOS, 
but it is at a distance to the K/K' points and therefore 
not captured by the effective theory. 

3. Armchair edge 

We now proceed with armchair type edges. According 
to Tab.lll the boundary projection matrix is given by 

Pa = ^ (1 -CT^ ® Ty) . (42) 

Then we obtain 

r.(fc) = 1 + ^{kBiy + ika,){l - P„) (43) 

and the surface Dirac-charge density reads 

/i.(fc,2/,) = 2P«Go(fc,-yO, (44) 

leading to [cf Eq. ^\ 

5Gi{k, Ui) = 2Go{k, y.^) iay P„ Go{k, -y.^) 

= -Go(fc,2/i)CT^Go(fc, -y,) ® Ty . (45) 

Surprisingly, since Ty is off-diagonal, the trace of 6Gi is 
zero and the boundary contribution to p in the armchair 
case vanishes. 

4- Infinite mass edge 

The calculation for the infinite mass edge is similar 
and for the surface Dirac-charge density we find as for 
the armchair case 

fi,ik,yi)^2P^Goik,-y,), (46) 

which leads to 

5Gi{k,yi) = ±Goik,y,)aMoik,~yi)®T, . (47) 

Similar as for the armchair edge, this expression is trace- 
less because Tr{Tz) ~ 0. However, we point out that 
even within individual valleys the boundary contribution 
to the DOS vanishes. This follows from the fact that 

oo 

J dy, Tr [Go(fc, y^)a,Go{k, -y,)] ^ (48) 



is an odd function of k and thus the corresponding in- 
tegral vanishes. This last fact has been already noticed 
by Berry and Mondragon^ for massless neutrinos in rel- 
ativistic billiards with infinite mass walls. 
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FIG. 3. (color online). Smooth part of the density of states 
for several graphene billiards with approximately the same 
area A ^ (140 a) ^, calculated numerically using a tight- 
binding code with only nearest neighbor coupling (solid lines). 
The numerical curves are obtained by first calculating exact 
eigenenergies and successive smoothing by replacing each en- 
ergy level by a Lorentzian with a half width at half maximum 
of 0.015 1. The dashed lines are the predictions of our the- 
ory, Eq. (gSll. From top to bottom: black: |9Vzz|/|(9V| = 1 
(zigzag triangle), blue: |9Vzz|/|9V| « 1/1.6 (Sinai shape), 
red: \dV^^\/\dV\ ~ 1/1.9 (rectangle), green: \dV^^\/\dV\ = 
(armchair triangle). 



D. Comparison with numerical results for various 
graphene billiards 

In summary, our result for the smooth DOS of a generic 
graphene billiard, neglecting the effect of next-nearest 
neighbors is 

p{k,)^-\k,\ + \dV,,\^S^{k,), (49) 

TT TT 

with I^V^^I being the total length of zigzag edges in the 
billiard. 

In Fig. 13] we compare our analytical result (|49p with 
results from numerical simulations for the graphene bil- 
liards shown as insets. For the numerical calculations we 
obtain the average DOS by computing eigenvalues of a 
corresponding tight-binding Hamiltonia n^^i^^ and subse- 
quent smoothing. All the billiards are chosen to have ap- 
proximately the same area. This is reflected in the com- 
mon slope of p for larger k^ , confirming the leading order 
term in the Weyl series. The different shapes and orien- 
tations give rise to different fractions of the zigzag bound- 
ary |i9Vzz|/|i9V|. While the boundaries of the equilateral 
triangles consist completely of either zigzag (black) or 
armchair (green) edges, both edge types are present in 
the rectangle (red) and in the non-integrable (modified) 
Sinai billiard (blue). We find very good agreement with 
our analytic prediction. We note that the dashed lines 
for the triangles and the rectangle do not involve any fit- 
ting, rather we have used the estimation /cmax = 7r/3a 
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FIG. 4. (color online). Smooth part of the density of states for 
the same systems as in Fig. [3] but with a relative next-nearest 
neighbor hopping strength t' =0.1. Solid lines show the nu- 
merical tight-binding results and dashed lines the predictions 
from Eq. (|50|) . For the smoothing we used Lorentzians with a 
half width at half maximum of 0.01 1. We used the same color 
coding as for Fig.O Inset: The tight-binding model exhibits a 
van Hove singularity at ks = — O.lt/fiVi? ~ —0.115 1/a. As a 
result the smoothed DOS shows a peak at the corresponding 
position (solid). 

from tight-binding theory. For the Sinai bilhard our the- 
ory allows to determine the total effective zigzag length 
= 516a. 

On the other hand, with nnn hopping we get from 
Eq. (1401) 

p-(fc.)^-|fc.| + |W..| ^~f^i^-^ (50) 
TT ZTrr 

In Fig.|4]we compare again this analytical result (dashed) 
with corresponding tight- binding calculations (solid). 
Also here we find good agreement with our analytic 
predicition for the surface term. Further towards the hole 
regime, i. e. to more negative energies, the tight-binding 
model has a van Hove singularity due to the edge state 
band edge at = — O.lt/ftwp ~ —0.115 1/a, as depicted 
in the inset of Fig.|4] (solid line). This peak is missing 
in our calculation, since in the effective Dirac theory the 
edge state dispersion is constantly linear for finite t' (cf. 
App. |B|). Note that also here, no additional fitting is 
involved (for the Sinai billiard we use |9Vzz| = 516 a ob- 
tained from the fit in Fig. [3]). 

From our discussion in this section it becomes clear 
that in principle the structure of a graphene fiake's 
boundary, i. e. the ratio between zigzag and armchair 
type edges, can be estimated from the behavior of the 
smoothed density of states at low energies. Hereby the 
formula predicts the spectral weight of the edge 
states dfefi pi(fci3) = |i9Vzz|/3a, which is model in- 
dependent, since the number of edge states is conserved. 
Note that Libisch tt al. have numerically investigated^! 
the averaged DOS of graphene billiards and found a p{kE) 



profile similar to that in Fig.|3l Related studies on edge 
states in graphene quantum dots have been performed in 
Ref. m 

IV. DENSITY OF STATES OSCILLATIONS 

A. The multiple reflection expansion in the 
semiclassical limit 

So far we have focused on the smooth part of the den- 
sity of states. In this section we study the oscillating 
part Pose- Our main result is an extension of Gutzwiller's 
trace formula^ to graphene systems with chaotic and 
regular classcial dynamics. We derive the trace formu- 
lae by evaluating Eq. ^ asymptotically in the semiclas- 
sical limit fc^L 3> 1. In other words we evaluate the 
boundary integrals in the MRE (fT6|) using the method 
of stationary phase. In the limit k^L 3> 1, the Han- 
kel functions become rapidly oscillating exponential func- 
tions of the boundary points. All other terms in Gn 
vary slowly along dV. Thus we evaluate them at the 
critical boundary points where the total phase of the ex- 
ponentials is stationary. There is another leading-order 
contribution to the boundary integrals that is of differ- 
ent origin, namely when the set of boundary points a = 
{a.N, ■ ■ ■ ,<y.i) leads to a singularity in the prefactor o^^-^^ . 
Due to the divergence of Go{a, f3) as |q; — /3| — > 0, quan- 
tum paths involving refiections at closely lying boundary 
points can give rise to such singularities. We show be- 
low that short range critical points occur only at zigzag 
edges. We treat these short range singularities at zigzag 
edges by resumming the MRE leading to a renormalized 
reflection operator. 

1. Resummation of short range processes 

The general method is outlined in Ref.[4^ Here we 
apply it to graphene. First we isolate the short range 
singularities: We define the action of an operator X on a 
function / 

i.f{a):= Jdaf3l{a,f3)f{f3). (51) 

dV 

In our case 

I(a,/3) = 2P„Go(a,/3)ia„^. (52) 

We now recast Eq. (|14p as 

fj.{a., x') = 2PaGo{a., x') — I/i(a, x') . (53) 

Furthermore we decompose T into a short range part Xj 
and a long range part T\: 

Is(a, /3) = X(a, 13) [1 - w{a - (3)] , 
Ix{oL,P)=I{oL,l3)w{a.- 13). ^ ' 
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Here ^(a — /3) is a smooth function, that is zero when- 
ever a is close to (3 and goes to one otherwise, so that 
integrating over f3 isolates the critical point (3 ~ a. This 
separation is however a formal one in that the specific 
form of w does not change the final result (see Ref.|6a for 
details). Then Eq. ^ leads to 

{l+is)^i{a,x') = 2PaGoia,x')^iifj.{a,x') (55) 

or with f = (i +Xs)"i 

fi{a, x') = 2fP„G'o(Q:, x') -fli ^{a, x') . (56) 

Now the renormalized Kernel Ii is free of short range 
singularities. Alternatively, in integral representation 

/x(a, x') = 2 J dap r(a, f3)PpGo{f3, x') (57) 

dV 

~ J dap J dap,T{a,l3)Ii{l3,l3')pi{l3',x'). 
dv av 

We note that the relevant structure of both terms in this 
expression is the same, since Ii contains the isolating 
function w and thus /3' can be considered to lie far away 
from /3 just as x' in the first term. In this way we have 
formally collected all the short range contributions in F 
and we are left with calculating 



2 / dapria,(3)PpGo{f3,x'). 



(58) 



dV 



We evaluate Eq. (j58|) again in the plane approximation 
and replace the boundary in the vicinity of a by a 
straight line in the direction of the tangent at a.. In 
our local coordinate system with x and y denoting coor- 
dinates in the tangential and normal directions, we ap- 
proximate a point f3 close to o; by /3 = {xp, yp) « (/3, 0), 
and write x' = (x',y') for a point x' far away from a 
(cf. Fig. [5]). Then the system is locally homogeneous 
along the straight boundary and we have 



(59) 
(60) 



F(a,/3)=F(a-/3), 

Go(/3,a;') = Go(/3-.T',-y')- 

In order to partial Fourier transform the expression (|58p , 
we use the convolution theorem to obtain (Pq, = Pp is 
constant along the straight boundary) 



d/3F(a-/3)P„Go(/?-a;',-y') 



^^g,fc(a-^')r(fc)P„Go(fc,-y')- (61) 



In fact we have calculated F(fc) already earlier, cf. 
Eq. (1201) and Eq. (gS]), leading to 



x' = {x',y') 




W,o) 

a = (a,0) 



FIG. 5. Notation in the local coordinate system spanned by 
the tangent and the normal to the boundary at a. Corrections 
to the approximation (3 ~ (/?, 0) are of subleading order in 
ksL, cf. IIIICll 



with the rcnormalizing factor 

( a(fc) 

1 



^,2 [a{k) ± kTz] for zz edges. 



i?«(fc) = <! ^ 

for ac and im edges . 

(63) 

We now define the renormalized free Green's function 
through its Fourier transform as 



Go(a,a;') 



dk 
2^ 



Mo^-^')ji^(^k)Go{k,-y'). (64) 



Finally we cast Eq. ([55)) for the charge layer fj, in position 
space into the form 

H{a, x') = 2F„Go(a, x') (65) 
- 2 dap PaGoioL,l3}w{a - f3) ia^i^ fi{(3, x'). 



dV 



The virtue of this equation is that it is free of short range 
singularities. 



2. Renormalized Green's function in the semiclassical limit 
With the definition 

k \ 



6{k) = arctan 



(66) 



we obtain from Eq. (|63p 

Ro.{k) -- 



cos[6{k)] e='=*f ('-•)'^^ for zz edges , 

1 for ac and im edges . 



(67) 

We compute Go (a, a;') in Eg. (|65p by performing the 
Fourier integral Eq. ([64)) [with from Eq. ([67)) ] within 
stationary phase approximation in the limit fc^L — > oo. 
We obtain the stationary phase point fco from 



F(fc)P„ = Rc.{k)Po, 



(62) 



dk 



kia-x')-VW^\y'\ 







(68) 
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yielding, in view of Eq. (|66[ 
tan[6'(fco) 



a ~ X 



\y'\ 



(69) 



The stationary phase point fco is such that the angle 0(fco) 
is equal to the angle that the vector x' — a includes with 
the normal at a, i.e. the classical angle of incidence. 
The stationary phase integration yields 



(70) 



Here Gf^ is the free Green's function in the semiclassical 
limit 



Gl%a,f3) 



2k, 



JkE\a~p\~iiT/A 



(1 + (7a,l3) 



41/ n\a- f3\ 

(71) 

where we use the short notation 
aa^p = (T ■ {a- j3)/\a- f3\ in Eq. ([71])- We note 
that expression ([71]) is closely related to the semiclassical 
Green's function for the free Schrodinger equation gg*^, 
namely 



Gl%a,f3) = k,gl%a,/3)il + aa.,p) 



(72) 



The matrix term reflects the chirality of the charge car- 
riers in graphene: the sublattice pseudospin is tied to 
the propagation direction and the projection (1 + a^/a) 
takes care of this. Eq. ((70|) together with Eq. ((65|) com- 
pletes our discussion of the short range divergencies and 
allows us to proceed with the long range contributions to 
the Green's function in the semiclassical limit. 



Note that the renormalization matrices Ra- account for 
possible short range singularities. 

Comparing Eq. with the MRE for the Hclmohz 
equation with Dirichlet boundary conditions^ shows 
that the scalar parts are very similar. The difference is 
that instead of factors ikEgQ^{cy.i^i, a.^), the MRE in Rcf. 
[39I has normal derivatives acting on the first argument 
cxi^i. In the semiclassical limit this leads to additional 
factors ifcjs cos(0i+i), where Oi^i denotes the angle be- 
tween the vector cc^+i — cti and the normal vector to the 
boundary at Q!i+i. We need not carry out the boundary 
integrals explicitly, but can immediately deduce 



G5^(a;,a;') = k^IxN Qwix, x') , 



where 



Kn = 



n^icos(0.) 



(76) 



(77) 



contains the pseudospin propagator as defined in 
Eq. ([73]), but a is now the vector of the classical re- 
flecti on p oints. The gj^{x,x') are well known, see e.g. 
Refs. [42I and IbtI. The stationary phase condition selects 
all sets of N stationary boundary points minimizing the 
phase aquired, and hence specifies classical trajectories 
of the system. We thus obtain our final expression for 
G^'^{x, x') in terms of a sum over classical trajectories 7 
that connect the points x' and x: 



G-(a;,cc') = — 



2~ 



7(a 



3. Semiclassical Green's function for graphene cavities 

In this section we evaluate the boundary integrals in 
the renormalized MRE in stationary phase approxima- 
tion. We consider the TV-reflcction term [cf. Eq. (|16p ] of 
the renormalized MRE, 

N 

GN{x,x')^{~2fY[ dao.,KN{a) fc^5g^(a;, ajv) 

. . .ikEgl''{a2,ai)ikEgo''{cxi,x') , (73) 

with a ~ (cci, ..CK;, ..ccjv)- In Eg. (j73|) we introduced 
the pseudospin propagator Ki^[a) that contains the 
graphene specific physics: 

N-l 

Kpfigt) = (1 + O-jj^ajv) Y\. '^ri^^Ra.Pa, (l + <7oLi + ^,a,) 
i=l 

xa„„^P„i (1 + T„^V)^(^) (74) 
with the separation function 

JV-l 



W{a) = Y[ w{ai+i - aj) . 



(75) 



Here, Lj, and are the length, the number of conju- 
gate points and the number of reflections at the boundary 
for the classical orbit 7. = K^., is the corresponding 
pseudospin propagator and 



dx± 
dp'j_ 



'1/2 



(79) 



measures the stability of the path 7 starting at x' with 
momentum p' and ending at x with momentum p. The _L 
denotes that the derivative involves only the projections 
perpendicular to the trajectory, which are scalars in two 
dimensions. 

Expression ([751) represents one main result of the 
present paper: The semiclassical charge dynamics for 
electrons and holes in a ballistic graphene flake is very 
similar to the case of electrons in Schrodinger billiards 
with Dirichlet boundary conditions. The graphene spe- 
cific physics is incorporated in the pseudospin dynamics 
described by K^. 

For a trajectory containing only one single reflection 
we have 



(1 + (^xa)a-n^'R'aPa{l + <y c 



(80) 
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Using the classical relations between the vectors x a 
and ex. ~ x' yields 

^ I e±*''^-crt„ (1 + (T^^,) for zz , ^^^^ 
1 e'^'^^CT^(l + o-Q.^,') foracandim. 

with 1/ according to Tab. HI With this result, we can ob- 
tain the pseudospin propagator for an arbitrary number 
of reflections by iteration. 



B. Trace formulae and semiclassical shell effects for 
classically integrable graphene billiards 

In this section we give two representative examples for 
trace formulae describing the oscillating part of the den- 
sity of states in graphene billiards that have classically 
integrable dynamics: circular and rectangular billiards 
with different types of graphene edges. We derive the 
corresponding semiclassical trace formula for the class of 
classically chaotic graphene cavities in Sec. IIV CI 

Orbits in regular systems are organized in families on 
classical invariant tori. An example of such a (periodic) 
orbit family is sketched for the circular billiard in Fig.|6l 
The members of a family possess the same classical prop- 
erties entering Eq. ((75)) such as action, length, stability, 
number of reflections and number of conjugate points. In 
order to compute the oscillatory part of the DOS from the 
semiclassical Green's function it is convenient to organize 
the trajectories in terms of tori, respectively families /, 
in the trace- integral, Eq. (|6]): 

p{kE) = -^JmJ^ JdxTr[Gf{x,x)] (82) 

/ Vf 

leading to the Berry- Tabor formula for posc hi terms of 
sums over families of periodic orbits organized on reso- 
nant tori^. The semiclassical pseudospin propagator for 
graphene does not alter the resonance condition (cf. the 
chaotic case IIV C[) , and for periodic classical orbits its 
trace Tr (K^) does not depend on the coordinates of the 
starting and end point: 

^ctix ^xot^ ^a^a^ ■ (^'^) 

Therefore, the integrals over V/ arc the same as for 
Schrodinger billiards with Dirichlet boundary conditions. 
Hence we can adapt the corresponding results by explic- 
itly including the correct pseudospin trace for each orbit 
family. 

The collective effect of orbit families giving rise to con- 
structive interference due to action degeneracies lead to 
pronounced signatures in the DOS of integrable systems 
known as shell effects^. We analyze below how such 
features are modified due to graphene edge effects. 




FIG. 6. Example of a family of degenerate classical orbits in a 
circular billiard. The black triangular orbit can be rotated by 
an arbitrary angle without changing its length. All resulting 
orbits contribute the same to the density of states, (adapted 
from Ref. [H.) 
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TABLE II. a) Energy levels knR of the circular billiard with 
infinite mass type edges obtained from the semiclassical trace 
formula Eq. (|90[l by summing over many classical orbits with 
^ = (TP) and by summing up all orbits approximately 
(TP (P)) Eq. H85|) compared to the quantum mechanical result 
(QM) Eq. (fSi)) . b), c) Energy levels k„L for square billiards 
with KL mod 2tt = 27r /3 {L — 200a "semiconducting") and 
KL mod 2tt = {L = 201a "metallic"), respectively. Again 
we compare the result from the semiclassical trace formula 
(|95p at 5 = with the quantum mechanical result (|C6I) . 

1. Circular billiard with infinite mass type edges 

We begin with a circular billiard with infinite 
mass type edges. Then the quantum energy levels 
Enm = ft-Vphnm are given by the intersections of Bessel 
functions^ 

Jn[kn,nR) = T J„+i (/s„„ii?) , (84) 

where R is the billiard radius, r = ±1 labels the two 
valleys and n^m & where m counts the intersections. 
For the semiclassical calculation of posc we adapt re- 
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suits for the Schrodingcr disk billiard as derived and dis- 
cussed in detail e.g. in Ref. [sgI . Periodic orbit families in 
the disk are labeled by the total number of reflections v 
and the winding number w, with v > 2w. Examples with 
w = 1, 2 are depicted in Fig.[7l We also allow for negative 
winding numbers w, and define the sign such that w > 
for clockwise going orbits and w < Q for anti-clockwise 
going orbits. Simple geometry gives for the length L^,w 
and the angle of rotation (p^.w aquired of an orbit (d, w) 



w 



Then the reflection angles read 



sgn(ti;) 



(85) 
(86) 



(87) 



Graphene physics enters through the pseudospin propa- 
gator, Eq. (|77p . with boundary matrix 

P« = (1 + r, ®c7t„)/2 (88) 

for the infinite mass case [see Eq. (|4]) and Tab.|T] ]. For 
an orbit {v,w) the trace over K yields 

TrA',,^ = i^Tr « a^e*"''"---) 

^ ^ /I \ I (—1)"/^ for even w , 
= AcosivOy nj) { (89) 
I for odd V . 

Equation ([SQ]) reveals the interesting property that only 
orbits with an even number of reflections are contributing 
to the oscillating DOS in the circular graphene billiard, 
while for odd v, the pseudospins are interfering destruc- 
tively. Note that this holds true also in each valley sep- 
arately, because in the case of odd v, the contributions 
from winding numbers w and —w have opposite signs. 

Adapting the expression for the circular Schrodinger 
billiard^i^ accordingly yields the semiclassical expres- 
sion for the oscillatory part of the DOS of the graphene 
disk: 



even 



(90) 



X sin3/2(^.„ „,)sin (^k^L„^^ + ^tt^ p-iiL,..u,/2f 

where = 1 if u = 2w and otherwise fy^w — 2. 

The last factor in Eq. (PH)) . giving rise to an exponen- 
tial suppression of orbits of length L^^^ > repre- 
sents a broadening of the peaks in the quantum density 
of states by convoluting p with a Gaussian of width ^. 
Such a broadening is additionally introduced to mimic 
e.g. temperature smearing or account for a finite life time 
of the quantum states, for instance due to residual dis- 
order scatteringiS. Thereby, Eq. ((90)) relates gross effects 
in smeared quantum spectra or experimental spectra ob- 
tained with limited resolution to the contributions from 
families of shortest periodic orbits^i^. 




to = 2 



FIG. 7. Classical periodic orbits representing families in the 
circular billiard, v is the total number of reflections along the 
orbit and w denotes the winding number. If {v,w) are not 
coprime the orbit is a repetition of a shorter primitive orbit. 
E. g. (4,2) is a repetition of (2,1) and (6,2) of (3,1). (Adapted 
from Ref. Igi.) 



Using the Poisson summation formula, we can approx- 
imately sum up the trace formula (j90p for ^ = and 
find the approximate eigenenergies kyw = xyw/R cor- 
responding to poles in the semiclassical sum, that fulfill 
the equation 

V + ^ = {2W + 1)[1 - arccos(VF/XyH^)/7r] 



2Xyw 



vw 



2W. 



(91) 



In Fig.[8]a)-c) we compare the results of the semiclas- 
sical trace formula (|90p with exact quantum results from 
Eq. ([84]) for the lower part of the graphene disk spec- 
trum. For ^ = [panel a)] even the exact quantum levels 
(blue circles) are reproduced with remarkable accuracy 
by the semiclassical theory [black peaks, see also numer- 
ical values in TabUJa)]. For every level, we have a sharp 
peak in the semiclassical result. An exception are the 
two levels close to fc^i? = 6, for which we have only one 
peak, though twice as high as the others, meaning that 
in the semiclassical expression the two levels are nearly 
degenerate. 

Panel b) shows the broadened spectrum for ^ = 0.3/i?. 
Again, the semiclassical result (solid line) is in very good 
agreement with the corresponding quantum result (dot- 
ted). For comparison, panel d) shows the same en- 
ergy range for the corresponding Schrodinger billiard. In 
Fig.[8]c) we have a closer look at which orbit families 
contribute. In fact we can see from Fig.[8]c) that the two 
shortest non-vanishing orbit families (2, 1) and (4, 1) al- 
ready yield a good approximation to the shell structure 
for e = 0.4/i?. 

Fig. ini shows the power spectrum of the exact quantum 
result (Gaussian convoluted with ^ = 0.4/i?). Evidently, 
only families with an even number of vertices v are con- 
tained in the spectrum, as semiclassically predicted. For 
example the triangular orbits (3, 1) that would give a 
peak at L/R = 5.2 and also the pentagram orbits (5,2) 
(i/i? = 9.5) do not contribute. The inset shows the same 
plot on a logarithmic scale, where the absence of the odd 
orbits is even more evident. 
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FIG. 8. (color online). Oscillating part posc of the density of 
states of a circular billiard as a function of ksR- a) Peaks are 
obtained from the semiclassical expression (|90p by summing 
up orbit families up to v, w = 400 for ^ = 0. Blue circles mark 
the positions of the exact quantum mechanical levels given by 
Eq. (|84p (See also Tab[TT]a)). b) Gaussian convoluted pose for 
^ — 0.3/ R. The full (dotted) curves show the semiclassical 
(quantum mechanical) results, c) Comparison between the 
full semiclassical orbit sum (dotted, ^ = 0.4/i?) with the con- 
tribution from the two shortest orbit families (2, 1) and (4, 1) 
(solid), d) Gorresponding results (for ^ = 0.4/i?) for a circular 
Schrodinger billiard with Dirichlet boundary conditions. 



2. Rectangular billiard with zigzag and armchair edges 

The rectangular billiard represents another promi- 
nent classically integrable geometry. While for the 
Schrodinger equation with Dirichlet boundary conditions 
this is a simple textbook problem, there is no explicit ex- 
pression for eigcncnergics of the graphene rectangle with 
two opposite zigzag and two opposite armchair edges. 
(For the derivation of a closed formula for the quantum 
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FIG. 9. Power spectrum of the Gaussian convoluted 
(^ = 0.3/ R) quantum density of states of the graphene disk 
with infinite mass edges. Peaks can be uniquely assigned to 
periodic orbit families {v,w), see text. Inset: Logarithmic 
respresentation. 




iN,M) = (0,1) 



FIG. 10. Families of periodic classical orbits in the rectangu- 
lar billiard. N{M) is the number of reflections at the bottom 
(left) side. 



cigencnergies in terms of a transcendental equation see 
App.[Cj. We will show that our semiclassical theory pro- 
vides a very good approximation to the quantum density 
of states. 

In the rectangle, the periodic orbit families can again 
be labeled with two indices. We denote by TV and M the 
number of reflections at the bottom zigzag (TV) and the 
left armchair (M) side of the rectangle with lengths 
and Ly respectively (see Fig.[TUl). The absolute values 
of the reflection angles at the zigzag and armchair edges 
then read 



l^z 



arctan 



ML, 
NL,, 



Ifilacl = I - = arctan (^^^ 



(92) 



From Eq. (|8ip we can read ofi' the following matrix factors 
for reflections with angles O^z and 0ac, respectively: 



•y 



(93) 



I a, lower zigzag edge, 

1 cTj, upper zigzag edge, 

"''^'^ left armchair edge, 

I a^e^^'^"'^' right armchair edge . 

This enables us to calculate the pseudospin trace of a 
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periodic orbit from family {N, M) as 

Tri^jVM = {-lf^cos{2MKL, - 2N%,\) . (94) 

This expression holds irrespective of the propagation di- 
rection along the orbit. Note also that the O^z in Eq. ([94|) 
occurs only due to the fact that we have different zigzag 
edges at the top and the bottom boundary (A- and B- 
terminated, respectively). Equation ([94| is now used 
to adapt the trace formula for the Schrodinger equation 
which has been derived e. g. in Refs. |5^ and l7ll to the 
case of graphcnc. Taking into account the interfering 
pseudospins in graphene, wc find 




NM 



vLnm 

,s (k^LNM - J) TrifjvM e-(«^"'^/2) 



(95) 



with length L^m = 'i^APLl + N'^L^ and TyKnm from 

Eq. dM]). Further /nm 1 if TV = or A/ = and 
otherwise /nm — 2. Note that the size of the bil- 
liard determines whether certain orbits contribute: The 
quantity KL^ can only take values that are multiples 
of 7r/3. In particular for KLx = mod 2tt—, families 
(TV, NLy/Lx) with odd N do not contribute according to 
Eq. (|94p . Further examples are the families (M, 0) and 
(0, N) for odd N and M respectively. They cancel each 
other exactly for KL^. = mod 2it because of the (—1)^ 
term in the pscudospin trace. 

In Fig. [TT] and TabHIlb). c) we compare the results from 



the semiclassical trace formula (|95|) for 



L 



with the quantum mechanical results obtained by solv- 
ing Eq. (|C6[) numerically. Again we find very good agree- 
ment with the quantum result. This is rather remarkable 
because of the complicated structure of the quantization 
condition (|C6p . The semiclassical predictions concerning 
the frequency content of the DOS oscillations are con- 
firmed in Fig.[TT] c) and d). For example the shortest 
orbits (1,0), (0,1) and (1,1) do not contribute for the 
system in d) [KL mod 2n = 0) due to destructive pscu- 
dospin interference, while they arc important in c) (KL 
mod 2n = 27r/3). 

Note that in TahM we find some additional levels 
from the semiclassical trace formula, which cannot be 
associated to quantum energy levels of the rectangle. 
Rather these peaks occur at positions that fulfill the 
quantization condition of a fictitious ID quantum well 
of width L with armchair boundary conditions. It is 
well known^ that this is an effect of sublcading order 
{[kEL]~-^/^ with respect to leading order) produced by 
orbits that 'graze' along the edges. 



C. Trace formula for classically chaotic graphene 
billiards 

Finally we consider classically chaotic graphene sys- 
tems. In this case no spatial symmetries are present that 



would give rise to an orbit degeneracy as in the regular 
case. From Eq. ([75)1 we already know that the final result 
differs from the trace formula for chaotic Schrodinger bil- 
liards only with respect to the pseudospin trace. Thus 
we have to work out how the spatial integral in Eq. (|6]) 
depends on this trace. To this end wc do not start di- 
rectly from the semiclassical Green's function ([75]) . but 
go one step back to Eq. ([75]) . In order to calculate the 
integral 



PJv(fcE) = — — 3^1 dxTr\GN(x,x)] 

TT J 



(96) 



we consider only the a;-dependent part of the integrand. 



In = I dx 



K{x,a) 



\/\x — a]y\\ai — x\ 



and choose the parametrization x = II + tt, where I is 
the direction from a^r to ai and t the direction perpen- 
dicular to I such that a right handed coordinate system 
results. The origin I = t = is at the point ccjv and we 
denote Ini = Iq^at — cn|. Then we can rewrite the phase 



ip{l,t)/kE = |a;-aAr| 



\ai — x\ 



Nl 



1 



2l\ 



'Nl 



(98) 



We are now evaluating the t-integral in stationary phase 
approximation assuming k^lNi 3> 1. The stationary 
phase point to is given by 



dip{l,to) k^lNito 



dt 1{Ini - 1) 
d^'f{l,to) kslNi 







dt^ 1{Ini-1)' 
ip{l, to) = ks Ini = kslaN - ai\ 



to = 0, (99) 

(100) 
(101) 



This means however that at the critical point to, the pscu- 
dospin propagator K{a) has no dependence on I left, 
since for t = Eq. (|83| holds. Thus the remaining inte- 
gral can be performed exactly: 



In = 



2TTl]\fl 



K{a}e 



iks let IV — ai 



(102) 



This tells us that as for the Green's function, we can 
essentially read off the result for p^^^ directly from 
the corresponding Dirichlet problem for the Schrodinger 
equation^ and find the Gutzwiller-typc trace formula for 
a chaotic graphene cavity 



PosAk.) = ^lHe^Tr(if,)A^e'^-^^ 



(103) 



Here the sum runs over all, infinitely many classical pe- 
riodic orbits 7, because the stationary phase points with 
t = to = are lying exactly on the straight line con- 
necting the last with the first reflection point, i. e. the 
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FIG. 11. Oscillating part of the density of states of a square billiard with two armchair edges and two zigzag edges {Lx — Ly — 
L). The left panels (a, c) show the results for a square with KL mod 2-k = 2n/3 ("semiconducting") and the right panels (b, 
d) are for a square with KL mod 27r = ("metallic"). Panels a) and b) show the Gaussian convoluted posc(fcB) for ^ = 0.3/L. 
The dotted curves represent the quantum mechanically exact results calculated with Eq. HC6P and broadened correspondingly. 
Panels c) and d) show the quantum mechanical power spectra (for ^ = 0.3). It is easy to identify the peaks associated with 
corresponding families {M,N). 



apperance of the pseudospin docs not affect the station- 
ary points. The classical amplitudes depend on the 
period, the stability and the Maslov index of the corre- 
sponding orbiii^. That means, except for h and the trace 
over Kj, accounting for the interference of pseudospins, 
the right-hand side of Eq. (|103p contains only classical 
quantities and has the same structure as Gutzwiller's 
trace formula. We note that in Ref. a semiclassi- 
cal trace formula is presented for posc, which however is 
not taking into account the boundaries required to obtain 
chaotic dynamics. Note that the expression ()103|) for posc 
is only valid for systems with isolated orbits, a prereq- 
uisite to evaluate the integral perpendicular to ocn ~ ai 
in stationary phase approximation. This is particularly 
fulfilled for chaotic systems. 

Expression (jl03[) allows, in principle, for computing 
semiclassical approximations for energy levels in chaotic 
graphcne billiards. Wc presume that this trace formula 
holds true more generally for classically chaotic graphene 
systems, not only billiards, with an appropriate gener- 
alization of the pseudospin evolution. Since the clas- 
sical dynamics of a graphene billiard is the same as 
that of a Schrodinger billiard, the convergence proper- 



tics of Eq. (|103p are expected to be similar to those of 
Gutzwiller's trace formula, with convergence problems 
linked to the exponential proliferation of periodic orbits 
with their length. In App. [DJ wc discuss the effect of 
weak bulk disorder on the trace formula (|103p . 

As Gutzwiller's trace formula for the case of quantum 
chaotic Schrodinger dynamics, the trace formula (|103p 
represents a suitable starting point to consider the sta- 
tistical properties of energy levels for chaotic graphene 
cavities, in particular universal spectral features within 
certain symmetry classes. Based on Eq. ()103|) we devote 
a major part of Refi^ to the semiclassical analysis of 
spectral statstics in graphene. There we will see that in- 
tervallcy scattering, semiclassically incorporated in the 
pseudospin dynamics, plays a key role for the effective 
symmetry class obeyed in graphene, e. g. unitary, or- 
thogonal or intermediate statistics between the two. 



V. CONCLUSION 

The growing ability to manufacture graphene-based 
nanostructures and their increasing role in the field of 
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graphene physics poses challenges to theory to treat con- 
finement effects. Addressing ballistic graphene cavities 
we have focussed on the effect of different types of edges, 
zigzag, armchair and inifinite mass type, on the spec- 
tral properties. The multiple reflection expansion used, 
combined with the semiclassical approximation, allows 
for incorporating and analyzing edge phenomena in a 
particularly transparent way, both for the mean density 
of states p as well as for the remaining oscillatory part: 
The leading-order Weyl contribution to p for graphene 
billiards scales with the phase space volume on the en- 
ergy shell, as for Schrodinger-type billiards. Edge ef- 
fects are expected to alter the perimeter correction to p, 
which is proportional to the total boundary length in the 
Schrodinger case with Dirichlet boundary conditions. We 
showed for graphene billiards that armchair and infinite 
mass edges do not give any perimeter contribution, while 
zigzag edges yield a characteristic low-energy term scal- 
ing with the length of the zigzag boundary. As analyzed 
in detail we could relate this boundary term in p to the 
average number of quantum zigzag edge states. Thereby, 
our approach allows for an alternative, analytical calcula- 
tion of the zigzag edge state contribution. For graphene 
nanostructures with unknown portion of zigzag-type edge 
segments, this enables one to estimate the effective zigzag 
edge length, respectively number of edge states, from the 
characteristic feature in p{E), see Figs. S] and [31 Hence, 
already the mean density of states of graphene flakes in- 
corporates important physical information. 

For the oscillatory contribution, Posc, to the density 
of states of graphene billiards we derive semiclassical 
trace formulae in terms of sums over classical periodic 
orbits. We show that, within the leading-order semi- 
classical approximation, the classical orbital dynamics 
entering into the semiclassical sums is the same as for 
Schrodinger billiards of the same geometry. This im- 
plies for regular graphene geometries Berry- Tabor like^ 
sums over families of orbits and for chaotic geometries a 
Gutzwiller type^ trace formula in terms of isolated pe- 
riodic trajectories. Edge effects enter into the contribu- 
tion of each periodic orbit (family) exclusively through 
the the pseudospin propagator and its trace along the 
orbit. This leads to a particularly transparent represen- 
tation of graphene edge phenomena. We gave a detailed 
interpretation for two representative regular systems: the 
graphene disk with inflnite mass edges and the graphene 
2d box with boundaries built from two zigzag and two 
armchair edges. The comparison with full quantum re- 
sults showed very good agreement, both for smeared 
spectra, highlighing the role of short, fundamental pe- 
riodic orbits, and on the level of individual energy levels, 
obtained semiclassically by summing up many orbit fam- 
ilies. 

A number of questions and further research directions 
is now arising from this work. They include the chal- 
lenge to generalize the semiclassical expressions for the 
density of states of clean billiards to cavities with im- 
purity scattering and systems with smooth conflnement 



potentials, more generally graphene with arbitrary clas- 
sical Hamiltonian dynamics, including also systems with 
mixed phase space. Second, the fact that our treat- 
ment of the zigzag edge associated average level den- 
sity proofs adequate for both settings, models without 
and with particle-hole breaking effects, e.g. from ncxt- 
nearcst-neighbor coupling, see Sec. IIII C 21 encourages to 
address zigzag edge magnetisroi^i^Zr— within this frame- 
work. Third, the semiclassical formalism developed al- 
lows for treating graphene nanostructures with bound- 
aries that can be viewed of being composed of many 
zigzag- and armchair-edge segments. In particular, ana- 
lytical expressions can be derived by treating long orbits 
with bounces off the different boundary segments in a 
statistical way. Fourth, the techniques used can be gen- 
eralized to quantum transport through open graphene 
nanostructures. 

In a second paper— we will particularly address the two 
last items and study spectral statistics (through the spec- 
tral form factor) of closed systems and transport prop- 
erties (weak localization, universal conductance fluctua- 
tions and shot noise) of open graphene billiards. 
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Appendix A: The singularity of a Dirac-charge layer 

Here we derive the expression (|13p inducting the dis- 
continuity of the Green's function at the boundarj*^. Us- 
ing the short distance asymptotic form for the Hankel 
function 



TT 



(Al) 



we obtain the short range singularities of the free Green's 
function from Eq. ^T0\\ 



f)(x,X > - 



i (T ■ {x — x') 



27r la; — x 



l\2 



(A2) 



If x' lies in the interior of V and cc is a point on the 
boundary 9V, 



lim Gq{x,x') — Go{a.,x') 



(A3) 
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is well defined and the first term in Eq. (|13p is trivially 
obtained from Eq. However if x' is on the boundary, 
the singular behavior of the Green's funetion becomes 
relevant. To see this, we perform the boundary integral 
in two parts, dividing dV into a small region Ds{cy.) = 
Cs{a)f] dV, where Cs{(y.) is a circle with radius S around 
a, and the remaining border -Da(a) = dV \ Dsia). We 
will take the limit (5 — > at the end of the calculation. 

We begin with the integration within Ds{a.). To this 
end we use the asymptotic expression for Gq and get 

^uwa) = lim lim dapGoix, f3)iani^^{f3,x') (A4) 

5—^0x^cy. J 



Ds{a) 



fTfi f (x 

— —u.(a.,x')cr- lim lim / daa: — 



-/3) 



where we took out of the integral and evaluated it at 
(3 = a.. Without loss of generality, we choose a — 0, 
X = \x\y and approximate Ds{cx) by a straight line along 
the X-axis, i.e. Ds{ot) = { f a; | f G [—(5,(5] }. Then we 
get 



a. 



>-Ds[a) 



— —n{oi, x')cr ■ lim lim / 



2tt 
(J. 



\x\y - 

-<5 



■—-^^[a.,x')cT-\\m lim 2 arctan((5/|a;|)y 

ZTT i5-s-0 |a;|->0 



1 



K'^,x') 



(A5) 



Since the kernel of the integral on Ds{oi.) has no singu- 
larity, it simply follows 

lim lim / dai3Go{x, l3)i(TngfJ,{(3,x') 

(5^0 X— s-ct J 



= J d(Ti3Go{a,f3)ian^fi{f3,x') . (A6) 

dV 

It is known from potential theory, that the integral on 
the right hand side exists^ and thus Eq. follows. 



Appendix B: Effective boundary condition for zigzag 
edges in the presence of next-nearest neighbor 
hopping 

It has been shown in Refs. [5^ and |6ll that the inclu- 
sion of next-nearest-neighbor (nnn) hopping in the tight- 
binding Hamiltonian of graphene has important conse- 
quences on the properties of the zigzag edge states. While 
for bulk graphene, up to a constant energy shift, the ef- 
fects are of subleading order in fc, for finite samples nnn 
hopping leads to an additional effective potential that 
is located solely on the edge atoms, therefore leading to 
qualitative changes of the edge state properties. These 
range from a finite dispersion to a complete change of the 
current profile in transporli^ . 



Here we neglect terms of higher order in k in the Hamil- 
tonian due to the nnn hopping and focus on the effects 
of the resulting edge potential. To this end we derive an 
effective boundary condition for the Dirac Hamiltonian 
with zigzag boundary. We consider a single zigzag edge, 
where the last row of atoms is located at i/q ~ a/\/3. Fur- 
thermore the graphene fiake shall be extended for y > yo, 
i.e. the last row of atoms is of B-type. The Hamiltonian 
is then given by^ 

t' 

H = Vj,(T p - hvp-5{y - ?;o)(l ~a^® t^) . (Bl) 

Here t' « 0.1 is the ratio of the next-nearest neighbor 
hopping constant, and the projection {l — az®Tz) ensures 
that the potential is located on the i?-sublattice. Similar 
edge potentials can model also adsorbants at graphene 
edges or edge magnctis m^^i^^ . 

The Dirac equation together with the Bloch theorem 
gives for the y-dependent part of the wavefunctions in 
the valley r = +1 



kEipsiy) = ki/jsiy) + 



dy 
dy 



(B2) 



-t'5{y-yo)ijB{y) 



(B3) 



Now we integrate these equations over a small window 
[?/o — yo + around the potential and take the limit 
e — )■ afterwards. Assuming that ip has at most a finite 
discontinuity at j/Qj we obtain from Eq. (|B2p 



(B4) 



lim iJBiy + e) - '4'B{yo - £) ^ , 



i. e. the _B— part of the spinor is continous. Thus we 
devide Eq. (|B3p by ipsiy) before integrating and get 



1 di'Aiy) 



yo+e 

= lim / - , ^ ^ 

£^0+ J ipB[y) oy 



lim 

£->0 + 



i'B{yo + £) ■fpB{yo-£) 



(B5) 



(B6) 



using integration by parts. For y < yo we employ the ac- 
tual zigzag boundary condition ipAiO) = 0, leading to the 
known expressions for the wavefunctions for y < yf^^ln^: 



-ipAiy) = Asm{qy) , 

V'B(y) = -T— [ikTsiniqy) + qcos{qy)] 



(B7) 



with longitudinal and transverse momenta k and q, re- 
spectively. Since the effective Dirac equation is valid for 
momenta that are much smaller than 1/a, we approxi- 
mate fc^a, qa, fca « to get 



lim 



ks sin(ga/\/3) 



0+ V'B(yo - e) ifc sin(ga/\/3) + 9Cos(ga/V3) 



«0, 
(B8) 
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which inserted into Eq. (|B6|) finally leads to the effective 
boundary condition 



IpA 



= t'. 



(B9) 



dV 



in agreement with a result found for similar edge poten- 
tials in Ref. [sil In an analogous way one can derive 
the effective boundary condition for the other valley as 
well as for A-terminated zigzag edges to end up with an 
effective boundary condition matrix 



it'ffy ± I'tz 



' <y.) (BIO) 



for all points at the edge ol. This expression turns into 
the usual zigzag matrix (|29p when t' = 0. 

We further derive the edge state dispersion and wave- 
function from the Dirac equation with the effective 
boundary condition 



V'A(o) = tVB(o). 



Be- 



(Bll) 
(B12) 



and 

xPB{v)^{Tk + dy)^A{y). (B13) 

For nonzero fc^ , Eq. (jBlip then leads to the condition 

A{kE - t'rk - it'q) = -B{kE - t'rk + it'q) . (B14) 

The bulk states result from this equation, when both 
sides arc nonzero. On the other hand the edge state 
results if this is not the case, e. g. fc^ — t'rk + it'q = 0. 
Solving this equation gives for negative rk the edge state 



,2fc 
» — 

k^ 



with the dispersion relation 

2kt'T 



k 



edge 



l+i'2 



2kt'T . 



(B15) 



(B16) 



This state exists only for negative (positive) momenta k 
in the valley K{K') (as for the case without nnn hopping) 
and has always a negative energy. 



Appendix C: Energy eigenvalues of a rectangular 
graphene flake 

Here wc present an implicit expression for the energy 
eigenvalues of a graphene rectangle with zigzag edges at 
y = Q and y = Ly and armchair edges at x = and 



X = Lx, respectively. To this end we start from a super- 
position of a forward and a backward propagating eigen- 
mode of an armchair nanoribbon with edges at a; = and 



■^{x,y)=A 



B 



( (g,„ - ifc)e'?"- \ 
kEe^i"-"" 
k^e-'i"-'' 

V {-qm + ifc)e-'"?-- J 

I {q^+lk)e'1^- \ 

V (-g,„ - ik)e~^i^- ) 



-iky 



(CI) 



where the qm are quantized according to 
mTT 



Lx 



— K m G Z . 



(C2) 



The spinors in Eq. (jCl|) arc solutions to the Dirac equa- 
tion when 



k' 



(C3) 



Now we impose the zigzag boundary conditions 
*A(a;,0) = *A'(a;,0) = -^B^x.Ly) = *B'(a;,Ly) = 0, 
which result in the two independent equations 

(g™ - ik)A + (q™ + tk)B = , (C4) 

g^fcL,^_^g-^fcL„^ (C5) 

These are solved for quantized knm that fulfill the tran- 
scendental equation 



qm tan(fc„m-^i/) 



(C6) 



With that we have formally solved the problem, the 
eigenenergies can be found e. g. by solving Eq. (|C6[) nu- 
merically. 



Appendix D: Effect of weak bulk disorder 

At this point wc briefly discuss the effect on the trace 
formula (|103p caused by smooth bulk disorder, which can 
be accounted for by an additional term 



(Dl) 



in the Hamiltonian, where V{x) is smooth on the scale 
of the lattice constant^^. In the semiclassical limit the 
Green's function for H + H' has been derived in Ref. 
without taking into account the boundaries. For the case 
of a Gaussian correlated disorder potential. 



{V{x)Vix')) = Coexp 



jx-x')"^ 
4A2 



(D2) 
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quantum calculations in the Boltzmann limit have been 
performed24iZ^. Under the assumption that the disor- 
der potential is weak enough that the classical trajecto- 
ries remain unaffected, we get for the impurity averaged 
Green's function^S 

{G'^^x,x')) « Gl''{x,x')cxp{^{5S^) /2h^) (D3) 

with 

X X 

{SS') = J^^;^ j dqj dq' {[{x - x') X VV{q)] 

x' x' 

y.[{x-x')y.W{q')]) /\x-x'\^ (D4) 
Kh^\x-x'\/l (D5) 

and the mean free path 

I - - . (D6) 



For smooth potentials the jump of the Green's function 
in Eq. (fT3| and hence also the MRE (|T6| remains un- 
changed, except that Go has to be replaced by its impu- 
rity averaged version (jD3p . Thus each summand in the 
semiclassical Green's function ((78)) for a graphene cav- 
ity aquires a damping factor e^^^/^'. In the trace inte- 
gral (|97| . these factors do not alter the stationary phase 
points, so that also in the trace formula in (|103p every pe- 
riodic orbit contribution is weighted with a factor e~^^/^' 
that improves convergence of the semiclassical trace for- 
mula. 
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